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^ . Abstract 

We investigate the relic density n x of non-relativistic long-lived or stable 
particles x m cosmological scenarios in which the temperature T is too low 
for x to achieve full chemical equilibrium. The case with a heavier particle 
decaying into \ is also investigated. We derive approximate solutions for n x (T) 
which accurately reproduce numerical results when full thermal equilibrium is 
not achieved. If full equilibrium is reached, our ansatz no longer reproduces 
the correct temperature dependence of the x number density. However, it does 
give the correct final relic density, to an accuracy of about 3% or better, for 
all cross sections and initial temperatures. 
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1 Introduction 



The production of massive, long-lived or stable relic particles x plays a crucial role 
in particle cosmology pQ. The perhaps most important example is the production 
of Massive Weakly-Interacting Particles (WIMPs), which may constitute most of 
the Dark Matter in the universe [2|. Alternatively, WIMPs may only be meta- 
stable, and decay into even more weakly interacting particles (e.g. gravitinos or 
axinos) that form the Dark Matter [3 j . Even if WIMP decays do not produce Dark 
Matter particles, the WIMP density is tightly constrained by analyses of Big Bang 
Nucleosynthesis (BBN)g]. 

It is usually assumed that the WIMPs were in full thermal and chemical equilib- 
rium in the radiation-dominated epoch after the period of last entropy production, 
which in standard cosmology means after the end of inflation. In this "standard" 
scenario the x number density n x (T) drops exponentially once the temperature T 
falls below the mass m x of the relic particles, until the freeze-out temperature Tp is 
reached, where the production of x particles from the thermal bath becomes negli- 
gible. In this case accurate semi-analytical expressions for n x (T <C Tp) have been 
derived one finds that the x re li c density is essentially inversely proportional 

to the thermal average of the effective % annihilation cross section into lighter parti- 
cles. The case of additional late entropy production, at T <C Tp, can also be treated 
analytically, by multiplying the standard result with a "dilution factor" due to the 
late-produced entropy [Zj. 

For typical WIMP scenarios, Tp ~ m x /20. The standard treatment can work 
only if the maximal temperature after inflation, usually called the reheat temperature 
Tr, is (much) larger than Tp. The assumption T R ^> Tp is not implausible, since 
the scale of inflation has to be quite high, typically ~ 10 13 GeV in simple models, in 
order to achieve the right order of magnitude of density perturbations jH]. On the 
other hand, we have direct observational evidence (from BBN) only for temperatures 
T < (few) MeV [HJEIIj which is well below Tp for most current WIMP candidates 
[2] . It is therefore legitimate to investigate scenarios with Tr < Tp [TTJ [T2J EH . 

We should emphasize at this point that Tr may not have been the highest tem- 
perature of the thermal plasma after inflation: given sufficiently fast thermalization, 
the inflaton decay products can attain a temperature T max 3> T R while the total 
energy density of the universe is still dominated by inflatons pQ. x particles may 
therefore have been in thermal equilibrium for some range of temperatures T > Tr 
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[TH |H1 E3 El EH! , even if they never were in equilibrium in the radiation-dominated 
epoch. However, an analytical treatment of the re-heating epoch where T > Tr was 
possible faces several complications not present in the radiation-dominated epoch: 
the entropy density was not constant, non-perturbative (and non-exponential) in- 
flaton decays might have been important , and there might have been significant 
non-thermal sources of x particles PB1E3E1- On the other hand, in supersymmet- 
ric scenarios thermalization of the inflaton decay products might be delayed by large 
vacuum expectation values of scalar fields along flat directions of the potential ^H] • 
In this paper we evade these complications by treating the x number density at some 
initial temperature To as a free parameter; in the absence of late entropy production, 
T should be close to the reheat temperature Tr (depending on the exact definition 
ofT fl ). 

Existing treatments of thermal WIMP production [5J El UH E3 E3 EEU EH] assume 
that n x had either achieved full equilibrium, or was completely out of equilibrium 
(i.e., annihilation of \ particles was always negligible). As already noted, in the 
former case one finds that the relic density is inversely proportional to the thermal 
average of the \ annihilation cross section. Not surprisingly, if \ annihilation can 
be neglected, one finds that the contribution to the \ relic density from thermal 
production is directly proportional to this cross section. Here we provide an approx- 
imate analytic treatment that also works in the intermediate region, where (for some 
range of temperatures) both thermal production and annihilation of x particles were 
important. It is based on an expansion in the effective annihilation cross section. To 
leading order, only the production term is kept in the Boltzmann equation describ- 
ing the evolution of n x (T); this corresponds to the "completely out of equilibrium" 
scenario analyzed previously. The first correction includes x annihilation, treating 
it as a small perturbation. This still allows an analytic solution, in terms of the 
exponential integral of first order Ei, which we only need for large values of its argu- 
ment. If n x (To) = 0, the first-order result is linear in the annihilation cross section 
a, while the correction is 0(a 3 ). Our most important, and (to us) rather surprising, 
result is that terms of higher order in a can be "re-summed" using a simple trick. 
This can be shown to be exact in the simple case where n x (To) > and thermal 
production of x particles is negligible*, and works numerically also for non-negligible 
thermal production. In fact, for T < T our formulae reproduce the exact numeri- 
cal results to 3% or better even for combinations of parameters where n x achieved 
*In this case the leading order result is trivial, i.e. O(a ), while the first correction is 0(a). 
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complete equilibrium, i.e. our new formulae are also accurate in scenarios where the 
"standard" result [3] is applicable. 

The outline of our paper is as follows. In Sec. 2 we briefly review the calculation 
of the relic abundance in the "standard" scenario, where it is assumed that the relic 
particles attained full thermal equilibrium. In Sec. 3 we will discuss our analytic 
calculation of the x re hc abundance in scenarios where the temperature was too low 
for x particles to have been in full equilibrium. In Sec. 4 we apply this method to 
more complicated scenarios, which include non-thermal x production from the decay 
of a heavier particle, still assuming the universe to be radiation dominated. Finally, 
Sec. 5 is devoted to a brief summary and some conclusions, while some technical 
details are given in the Appendix. 



2 Relic Abundance in the Standard Cosmological 
Scenario 

We briefly review the calculation of the relic density of long-lived or stable particles 
X in the standard cosmological scenario [JJJ, which assumes that the relic particles 
were in thermal equilibrium in the early universe and decoupled when they were 
non-relativistic. The relic density can be calculated by solving the Boltzmann equa- 
tion which describes the time evolution of the number density n x in the expanding 
universe [T], 

dn 



^ + 3ffn x = -H( n ;-n;j , (1) 

with n x eq being the equilibrium number density of the relic particles, H the Hubble 
parameter and (av) the thermal average of the annihilation cross section a multiplied 
with the relative velocity v of the two annihilating x particles. The first (second) 
term on the right-hand side (rhs) of Eq.((TJ) describes the decrease (increase) of the 
number density due to annihilation into (production from) lighter particles. The 
equilibrium density in the non-relativistic limit is given by 

n x , eq = 9 x {^P)' 2 ^' T > ( 2 ) 

where m x and g x are the mass and the number of internal degrees of freedom of 
X, respectively. In the standard cosmological scenario, it is assumed that x was 
in thermal equilibrium for T > m x . In other words, x rapidly annihilated with its 
own antiparticle into lighter states and vice versa. At later times T <C m x , the 
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annihilation rate T x = n x (av) dropped below the expansion rate H. Therefore x 
particles were no longer able to annihilate efficiently and the number density per co- 
moving volume became constant. The temperature at which the particle decouples 
from the thermal bath is called freeze-out temperature Tp. 

The Boltzmann equation (JJJ can be rewritten by introducing the new variables 
Y x = n x /s and K^eq = n Xfiq /s, where the entropy density s = (27r 2 /45)g*T 3 with g* 
being the number of the relativistic degrees of freedom. Assuming that the universe 
expands adiabatically, the entropy per comoving volume is conserved. Hence we 
obtain^ n x + 3Hn x = sY x . In the radiation dominated era the Hubble parameter is 
given by 

H=*HJ», t=-L, (3) 
Mpi V 90 ' 2H ' v ' 

where Mpi is the reduced Planck mass, Mpi = 2.4 x 10 18 GeV. By introducing the 

inverse scaled temperature x = m/T, the Boltzmann equation (JTJ) becomes 

dY 

= -1.32 m x M vl ^gZ{av)x- 2 {Y^ - Y 2 J . (4) 

In most (although not all cases the cross section is well approximated by a 
non-relativistic expansion: 

(av) = a + b(v 2 ) + 0({v 4 )) = a + 6b/x + 0(l/x 2 ) . (5) 

Here a is the v — > limit of the contribution to av where the two annihilating % 
particles are in an 5" wave. If S wave annihilation is suppressed, b describes the P 
wave contribution to av. In the following we treat a and b as free parameters. In 
terms of the variable A = Y x — Y x eq , the Boltzmann equation (@J) can be rewritten 

as 

where 

A = 1.32 m x Mpiy/g^(a + 6b/x) x~ 2 . (7) 

An analytic solution can be obtained by considering the equation in two extreme 
regimes. At early times (x -C xf), Y tracks its equilibrium value Y eq very closely. 
Therefore A and dA/dx are small. Ignoring A 2 and dA/dx, we obtain 

A ' 2.\ * ^ (8) 

^Here we assume — 0. This is usually justified since, as we will see below, n x has non-trivial 
time dependence only for a rather narrow range of temperatures; moreover, except during the QCD 
phase transition at T ~ 200 MeV, g* changes slowly, i.e. dg*/dx -C <?*• 
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where we used dY XtEq /dx ~ — Y Xfiq for x ^> 1. At late times (x ^> xp), one can 
ignore the production term in the Boltzmann equation: 

dA , * 9 

— ~ -AA . 9 

ax 

Integrating this equation from to infinity and using the fact that A(xf) 3> A(oo), 
we have 

y Xi00 = y x (x > xjr) = XF . (10) 

1.32 m x MpiA/g*(x_F)(a + 
It is useful to express the energy density as f2 x = p x j p c , where p c = 3HqMp { is the 
critical density of the universe. The present energy density of the relic particle is 
given by p x = m x n Xt00 = m x soY XjOC , with so ~ 2900 cm -3 being the present entropy 
density. Finally, we obtain the standard approximate formula for the relic density: 

, _, 8.7 x 10-" *, GeV^ (n) 
\/g*(xp)(a + 3b/x F ) 

where h is the scaled Hubble constant, h ~ 0.7. Notice that the relic density of the 
particle is inversely proportional to the annihilation cross section and that there is 
no explicit dependence on the mass of the particle. Calculating the cross section and 
the freeze-out temperature is sufficient for predicting the relic density. Freeze-out 
occurs when the deviation A is of the same order as the equilibrium value: 

A(x F ) = OWa*) , (12) 

where £ is a numerical constant of order unity. Substituting the early time solution 
of Eq.(jHJ) into this equation, xp is obtained by iteratively solving 

0.382 £m x M Pl g x (a + 6b/x F ) 

Xp — in = . [L6) 

x/xpg^xp) 

It is known that the choice £ = v2 — 1 gives a good approximation of exact nu- 
merical results for the relic density (|llj). The decoupling temperature depends only 
logarithmically on the cross section. For WIMPs, we typically obtain ~ 22. 



3 Relic Abundance in a Low— Temperature Sce- 
nario 

Eq. (jllj) implies that the relic density predicted in the standard cosmological scenario, 
in which x particles are assumed to have been in full equilibrium, would be quite 
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high unless the cross section is as large as* ~ 10~ 9 GeV -2 . Bearing this situation in 
mind, it is important to explore scenarios where the relic density comes out smaller 
than the standard calculation and find a useful formula which properly describes the 
behavior of the relic abundance. 

For later convenience we first rewrite the Boltzmann equation (J3J), using Eq.(j2|: 

w = -'{ a+ T)h<X-^' (14) 

where 

/ = 1.32V£m x Mpi, c = 0.0210 (15) 

are constants. Eqs. (@J) and ()14j) assume that x remains in kinetic equilibrium through 
the entire period with non-negligible time dependence of Y x . This is reasonable, 
since kinetic equilibrium can be maintained through elastic scattering of \ particles 
on particles in the thermal plasma. The rate for such reactions exceeds the \ annihi- 
lation rate by a factor oc Y~ x > 10 7 for temperatures of interest. For our numerical 
examples, we consider a Majorana fermion with m x = 100 GeV and g x = 2 as the 
relic particle. We choose the relativistic degrees of freedom to be g* = 90; this ap- 
proximates the prediction of the Standard Model of particle physics for temperatures 
around 10 GeV. 

Figure ^ shows that the relic density can be reduced if the particles never reach 
thermal equilibrium because of the low reheat temperature after inflation. The solid 
red curves depict the predicted present relic density fl x h 2 as function of a (a) and b 
(b) defined in Eq.fjSJ). Here we assume that the relic abundance vanished at the initial 
temperature of xq = 22, which is around the typical WIMP decoupling temperature. 
Here, as well as in the subsequent figures, the "exact" numerical solution of the 
Boltzmann equation (|14j) has been obtained using the Runge-Kutta algorithm, with 
a step size that increases quickly with increasing x — xq. For large cross section we 
observe Q x h 2 oc l/(av), in accord with the "standard" prediction (|llj). However, 
when the cross section is reduced, the relic density reaches a maximum, and then 
decreases oc (crv). For the given choice of initial conditions, there are therefore two 
distinct ranges in (crv) where the relic density comes out in the desired range |20j . 

In the following we attempt to find a convenient analytic formula applicable even 

to low temperature scenarios. As zeroth order solution of Eq. (jl4|) we consider the 

*We use natural units, where H = c = feg = 1, so that both a and av have dimensions GeV -2 . 
Numerically, 10 -9 GeV" 2 = 0.388 pb = 1.16 • 1(T 26 cm 3 /s. 
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Figure 1: Predicted present relic density £l x h 2 as function of the a and b contributions to the 
total cross section, see Eq.©; in frame (a), b — whereas in (b), a ~ 0. We consider two extreme 
cases: \ particles were in full thermal equilibrium (dotted blue line) or the number density of x 
vanished (solid red line) at Xq = 22. The two horizontal double-dotted black lines correspond to 
the la upper and lower bounds of the dark matter abundance |2U|. 



case where \ annihilation is completely negligible 

dY 



dx 



fc(ax + 66) e~ 



2x 



This equation can easily be integrated, giving 



Yo(x) = fc [-(x e- 2xo - 



For x ^> xq, the relic abundance of the particles becomes constant, 



Yo,oo =Y (X^> Xq) = fc 



-x e 



-2x 



1^36)' 



The corresponding prediction for the present relic density is given by 

tt x h 2 = 2.8 x 10 8 m x F 0jOO GeV" 1 . 



(16) 



(17) 



(19) 



Notice that the relic density is proportional to the cross section, although the coef- 
ficient of proportionality depends on whether a or b is dominant. 

So far no analytic solution has been known for the in-between case where both 
annihilation and production play a crucial role in determining the relic abundance 
while thermal equilibrium is not fully achieved. We now attempt to connect the 
standard scenario (T R > T F ) and the low reheat temperature scenario (T R < T F ) 
using some analytic method. 
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Since we already have the solution only including the production term, the most 
natural extension is to add a correction term which describes the effect of annihilation 
on the solution for the pure production case: 



Y 1 = Y + 5. 



(20) 



By definition 5 vanishes at the initial temperature. Since it describes the effect of x 
annihilation, it is negative for x > Xq. As long as \S\ is small compared to Y Q , the 
evolution equation for 5 is given by 



d8 
dx 



-f a 



66 \ Yn(x) 2 



xl or 



(21) 



Using Eq. ffTTj) for Y (x), this can again be integrated: 
S(x) = -fc 2 



^a 3 F 4 (x, x ) + ^a 2 (a + 18b)Ff(x, x ) 

+^a(a + 126) (a + 366)F 2 4 (x, x ) + -b(a + 126) 2 F 3 4 (x, x ) 
lb o 



a 2 F 2 (x, x Q ) + -a(a + 24b)F 2 (x, x ) + 36(a + 126)F 3 2 ( 



+ ^0,oo/ 2 C 

- Y 2 t00 f [aF 2 °(x, x ) + 66F 3 (x, x 



where 



F™(x,x ) 



x g— mi 

dt — — , m = 0,2,4, n = l,2,3 



t' 



(22) 



(23) 



The functions F™(x, Xq) can be expressed analytically in terms of the exponential in- 
tegral of first order Ei(x); a complete list of the relevant F™ is given in the Appendix, 
Eqs.(JH2). At late times, x — > oo, this simplifies to 



S(x — > oo) 



3_2_-4x 



-/Ve 



a 2 (a + 606) 9ab(a - 166) 



T x ° + 16 
96(5a 2 - 56a6 + 966 2 ) 



ix 



- f 2 ce- 2 *°Y x (x ) 



32x 2 

' 2 9a6 96(a - 46) 



•i'o 



2*2 



(24) 



\Xq Xq 

where we omit higher order terms than 0{1/xq). Notice that we discard Oi\jx 2 
and 0(l/x 3 ) terms in (av), which also contribute to higher order terms in Eq. 
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If a ^ we therefore expect additional terms 0(1/ xq) from terms not included 
in Eq.flHJ); if a = 0, higher order terms in the expansion of the cross section only 
contribute at 0(l/xl) in Eq.fl2U). 

Since, for vanishing initial abundance, Y is proportional to the cross section a, 5 
is proportional to a 3 . On the other hand, for sufficiently large cross section we want 
to recover the standard expression, where Y x (x — > oo) oc l/(av). This suggests to 
rewrite our ansatz ([211)) as 

Y >= Y ° +s = Y °( 1+ v^T^wr Y ^ < 25) 

Although the final approximate equality in Eq.(j2SJ) only holds for \5\ <C Y , we 
note that the resulting expression has the right behavior, Yx >r ocl/a, for large cross 
section. In the following we will show that this "resummation" of the correction 5 
is indeed able to describe the relic density for a wide range of cross sections and 
temperatures, including scenarios where the standard treatment is applicable. 

In fact, this ansatz solves the Boltzmann equation (fTlj) exactly in the simple 
case where thermal x production can be ignored, but Y x (x ) is sizable, leading to 
significant \ annihilation. In this case Eq. (jl4|) reduces to 



— * = (a+- -J. 26 

ax \ x ) x A 

This equation can easily be solved analytically. The solution decreases monotonically 
from its initial value Y x (xq): 

Y = y *(^o) (27) 

x l + fY x (x )[a(l/x -l/x) + 3b(l/x 2 -l/x 2 )]- 1 } 

In order to treat this case using the formalism of Eqs. (jl6|) - (}2l)j) . we simply drop all 
terms which depend exponentially on x or xq\ these terms come from thermal x 
production, and are obviously very small for sufficiently small initial temperature. 
The zeroth order solution (fTTj) then obviously reduces to the constant Y x (xo), and 
the correction S of Eq. (j22J) simplifies to 



(28) 



8(x) -> -f(Y x (xo)) 2 [aF°(x,xo) + GbF°(x,x )] 



' /l 1\ / 1 1 
a + 36 \ —z 

\Xq X J \Xq X z 



in the last step we have used the last two Eqs. pfij) . Inserting this in the last expres- 
sion in Eq.(|23j). we indeed recover the exact solution (|2Tjl . as advertised. 
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In principle, we can add further correction terms to the first order approximation 
of Eq.flHJ), 



The above discussion shows that this corresponds to an expansion in powers of (av) . 
Since Yq > and 5 < by definition, the systematic expansion will lead to an 
alternating series which possesses good convergence properties. However, this type 
of expansion is quite cumbersome because \S\ often dominates over Yq for not very 
small cross sections, as we will explicitly see later. Therefore the re-summed ansatz 
Y\, r of Eq. (j25)l is much more convenient. We will see that it often provides a good 
approximation to the exact solution even if thermal \ production is not negligible. 

In Fig. El we present the evolution of the exact, numerical solution Y x (solid 
red), Y 1>r (dotted blue), Y Xfiq (double-dotted black) and \5\ (short-dashed violet) as 
function of x — xq. Here we consider vanishing initial \ density, Y x (xq = 22) = 0. 
Clearly the first order approximation Y\ of Eq. (}2T)|) fails to reproduce the exact result 
once \8\ becomes comparable to Yq. On the contrary, frames (a) and (c) show that 
the re-summed ansatz Y l r of Eq. (J25j) reproduces the numerical solution very well 
for all x > Xq if a < 10~ 9 GeV -2 and b < 10~ 8 GeV~ 2 . However, for intermediate 
values of x — Xq, the disagreement between Y^ r and the exact solution becomes large 
as the cross section increases. In frames (b) and (d) of Fig. El sizable deviations from 
the exact value are observed at x — Xq ~ 1 for a = 10~ 8 GeV -2 or b = 10~ 7 GeV~ 2 . 
For larger x the deviation becomes smaller again, and for x 3> Xq the difference is 
insignificant even for these large cross sections. 

We also analyzed scenarios with sizable initial x abundance, Y x (xo) ^ 0. Figure El 
shows that the re-summed ansatz again matches the numerical result very well for 
all values of x if a < 10 -9 GeV~ 2 . This is not surprising since, as we saw in the 
discussion of Eq.(|2S|l. it reproduces the exact solution if Y x (xo) dominates over the 
thermal contribution. For a = 10~ 8 GeV -2 , Yi j7 . again starts to deviate from the 
exact numerical solution at x ~ 0.1, but approaches it for x ^> xq. Note also that 
already for the smaller cross section chosen in this Figure, the final relic density is 
almost independent of Y x (xq). 

Let us take a closer look at the difference between the exact solution and the 
re-summed ansatz. To this end, we define the deviation e by 



Y x = Y + 5 + 5 2 + 8 3 + ■ ■ ■ . 



(29) 



1 - S/Y 



+ e. 



(30) 
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(c) (d) 

Figure 2: Evolution of the exact solution Y x (solid red curves), Yt <r of Ea. (|2"5)l (dotted blue), the 
equilibrium density Y x ^ cq of Eq.JSJ) (double-dotted black), and |<5| of Ea. l22fl (short-dashed violet) 
as function of x — Xq. The initial abundance is assumed to be Y x (xq = 22) = 0. We take (a) 
a = 1CT 9 GeV~ 2 , b = 0, (b) a = 1CT 8 GeV" 2 , b = 0, (c) a = 0, b = 1CT 8 GeV" 2 , and (d) a = 0, 
b = 10~ 7 GcV~ 2 . In frames (a) and (c) the curves for Y\, r practically coincide with the solid lines. 



Inserting this ansatz into the Boltzmann equation (jlj) leads to the evolution equation 
for e: 

dx x 2 

which again resembles the Boltzmann equation. Since initially e = 0, our re-summed 
ansatz works very well as long as 5/Y remains suppressed. Note that the inhomo- 
geneous term on the rhs of Eq.([3ip is of order (5/Y ) 2 . The analogous correction 
to our original first order solution Yi of Eq. (|2T)|) would start at O(5/Y ). Since 
this inhomogeneous term is positive, e(x) > for all x > xq, i.e. Yi >r , like Yi, al- 
ways under-estimates the exact solution. As |<5|/Yq grows, the last term in Eq.(J3T} 



e 2 + 2e 



Yp (s/Y y 

- S/Y (1 - 5/Y,y 



(31) 
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(c) (d) 

Figure 3: Evolution of Y x (solid red curves), Y\. r (dotted blue), Y x , oq (double-dotted black) and 
\S\ (short-dashed violet) as function of x — xq. Here we take (a) a — 10~ 9 GeV" 2 , Y x (xo) — 10 -8 , 
(b) a = 1(T 9 GeV- 2 , Y x (x ) = 1(T 10 , (c) a = 1(T 8 GoV" 2 , Y x (x ) = 10~ 7 and (d) a = 1(T 8 
GeV -2 , Y x (xq) — 10 -10 . The other parameters are as in Fig.El 



can become sizable. Note, however, that it is multiplied with (Y Xteq ) , which drops 
oc exp (— 2x) with increasing x. Therefore e becomes large only if \S\ reaches values 
of order of Y for x — Xq < 1 . The homogeneous terms in Eq. (p?Tj) imply that for large 
x — Xq the deviation e decreases again, similar to the WIMP relic abundance Y x . This 
situation is depicted in Fig. EJ which shows the evolutions of (upper curves) 

and e/Y x (lower curves) as function of x — xq for a = 3 x 10~ 8 GeV~ 2 (solid red), 
a = 1(T 8 GeV" 2 (dotted blue) and a = 3 x KT 9 GeV" 2 (double-dotted black). Here 
we choose 6 = and Y x (xq = 22) = 0. Even in the case where e becomes sizable for 
intermediate values of x, it eventually diminishes and hence our analytical formula 
succeeds in reproducing the present relic abundance Y x (x — ► oo) fairly well. 
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x-x n 

Figure 4: Evolution of |<5|/Y^ (upper curves) and e/Y x (lower curves) as function of x ~ xq 
for a = 3 x 1CT 8 GeV" 2 (solid red), a = 1CT 8 GeV" 2 (dotted blue) and a = 3 x 1(T 9 GeV~ 2 
(double-dotted black). Here we choose 6 = and Y x (xq = 22) = 0. 

Let us turn to a discussion of the dependence of the present relic abundance 
on the initial temperature. In Fig. we plot the present relic density evaluated 
numerically (solid red curves), the old standard approximation (dotted blue) and 
our new approximation (double-dotted black) as function of xq. Here we take (a) 
a = 1(T 8 GeV~ 2 , b = and (b) a = KT 9 GeV~ 2 , 6 = 0. We find that our 
approximation agrees with the exact result very well for x > x f- On the other 
hand, for xq < Xp, our approximation gives too small an abundance^ while the old 
approximation works very well. The transition between the two regimes is very 
sharp. For xq = xf + 2, the old approximation over-estimates the relic abundance 
by as much as an order of magnitude, while for Xq — X p both the old and the new 
approximation work well. 

We found that for vanishing initial x density, Y x (xq) = 0, different values of the 
cross section lead to a universal behavior when the present relic density is expressed 
as function of Xq — xp and in units of the relic density for xq <C xp. This can be 
seen from the analytic solution we have obtained. For xq xp it is obvious that 
Q x (xo)/Q x (xo xf) is nothing but unity and independent of the cross section. For 
Xq ^> Xf, the exact solution is roughly given by the zeroth order approximation Yq, 
which scales like ae~ 2x ° if a dominates and the initial abundance vanishes. Mean- 

TFor xq <§; xp, our expressions predict fl x h 2 oc Xq- 
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10 15 20 25 30 10 15 20 25 30 

x x 

(a) (b) 

Figure 5: The present relic density evaluated numerically (solid red curves), the old standard 
approximation (dotted blue) and our new approximation (double-dotted black) as function of xo- 
Here we take (a) a = 1CT 8 GeV" 2 , 6 = and (b) a = 1CT 9 GcV~ 2 , 6 = 0. 




Figure 6: n x (x )/£l x (x <C xp) as function of x ~ x F . In the left frame, a = 10~ 7 GeV -2 
(solid red curves), 1CT 8 GeV" 2 (dotted blue) and 1(T 9 GcV~ 2 (double-dotted black) with 6 = 0, 
whereas in the right frame, 6 = 10~ 6 GeV~ 2 (solid red), 10~ 7 GeV~ 2 (dotted green) and 10 -8 
GeV -2 (double-dotted black) with a = 0. 



while, Eq. ()13|) shows that Xp is roughly proportional to In a. Therefore we obtain 
the relation 

^ (Xo) oc ^ oc e- 2 (— ) , (32) 



f2 x (x -C x F ) 1/a 

which has no explicit dependence on the cross section. The same argument is ap- 
plicable to the case where b is dominant. In Fig. |U]we plot the ratio of the exact 
present relic density to the value for x <S xp, Q x (xq)/Q x (xo <C xf), as function 
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of xq — xf for various values of a and b. These figures clearly show the expected 
scaling behavior both for a ^ 0, b = (left frame) and for a = 0, b ^ (right frame). 
However, for Y x (xq) ^ 0, no such scaling exists, apart from the fairly obvious result 
that Y x (x ^> Xq) becomes independent of Y x (x ) if x <C xp. 

Fig. El shows that Yi )T (xq, x — > oo) has a well defined maximum when xq is 
varied. This maximum occurs at a value £o,max which is close, but not identical, to 
the decoupling temperature xf of Eq. (jl3|) . From the asymptotic expressions for Yq, 
Eq.fEEJ), and 5, Eq. pijl . we find for Y x (xq) = 0: 



•^O.max — „ In 



1 / 2 c(a + 66/x , 



\2 



max ^ 



max 



2 4x 0: 
_ ln 0.096 m x M Pl g x (a + 66/x 0;max ) 

In deriving this equation, we neglect non-leading terms in l/xo,max in each combi- 
nation of a and b} Notice that x , m ax coincides with xp of Eq. (fT3|) . if one chooses 
£ = 1/4 (rather than £ = \/2- 1). 

Since the actual relic density is already practically independent of x for x < 
%o,max we can construct a new semi-analytic solution which describes the relic density 
for the whole range of x$: for x > ^o,max, compute the relic density from Yi tr (xo), 
but for xq < x 0) max, use Y l r (x 0) max) instead. 

The ratio of this semi-analytic result Q\ tr to the exact value Q x is depicted in 

Fig. [7| As noted earlier, our approximation becomes exact for xq > xp. For smaller 

Xq the new approximation still slightly under-estimates the correct answer, but the 

deviation is at most 1.7% for 6 = (left frame), and 3.0% for a = (right frame). 

On the other hand, in the same region the old standard approximation reproduces 

the present relic abundance within 1% error. We thus see that for xq < xp, this 

new expression works nearly as well as the old standard result^ of course, the old 

result fails badly for xq > xf- Finally, since by definition Yi r depends only weakly 

on xq for xq ~ xo, m ax, the latter quantity need not be calculated very precisely; in 

practice, setting Xo, m ax — 20 in the rhs of Eq. (J33|) is often sufficient. In contrast, the 

■'"The next-to-leading correction to the pure a— term would have been relevant, but it cancels. 
The non-leading corrections to terms that require both a and b to be non-zero are numerically 
insignificant, and of the same order as terms omitted in the expansion JSJ) of the annihilation cross 
section. 

^However, if a = 0, we should expect O(10%) corrections to the relic density from higher order 
terms in the expansion © of the cross section; if a ^ 0, these higher order terms should only 
contribute 0(1%). 
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(a) (b) 

Figure 7: Ratios of approximate and exact results for the relic density f2i r /f2 x as function of 
— ^o.max, for a ^ 0, b — (left frame) and a = 0,6 ^ (right frame). The curves use Y\ iT 
with xq replaced by max(xo, £o,max), see Ea. (|55|) . In the left frame, a — 10~ 8 GeV~ 2 (solid red 
curves), 1CT 9 GeV~ 2 (dotted blue), 1(T 10 GeV~ 2 (double-dotted black), 1CT 11 GeV~ 2 (short- 
dashed violet) with 6 = 0, whereas in the right frame, b = 10~ 7 GeV -2 (solid red), 10 -8 GeV~ 2 
(dotted blue), 10~ 9 GeV" 2 (double-dotted black), 10~ 10 GeV" 2 (short-dashed violet) with a = 0. 



standard approximation (|TT|) depends linearly (for b = 0) or even quadratically (for 
a = 0) on xp; several iterations are therefore required to solve Eq. (fT3j) to sufficient 
accuracy. Altogether, our new semi-analytic formula is evidently a quite powerful 
tool in calculating the density of cold relics. 



4 Relic Abundance Including the Decay of Heav- 
ier Particles 

In this section we investigate a scenario where unstable heavy particles decay into 
long-lived or stable particles x- We assume that decays out of thermal equilibrium, 
so that production is negligible; however, we include both thermal and non-thermal 
production of x particles. For example in some supersymmetric models neutralinos, 
which are stable due to R-parity, can be produced non-thermally through the decay 
of moduli [21] or gravitinos after the end of inflation. The number densities of x an d 
obey the following coupled Boltzmann equations: 
din 

-^ + 3Hn x = -(av)(n 2 x -nl cq ) + NT^n^ 

^± + 3Hn ^ = _ r ^, (34) 
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where N is the average number of x particles produced in a decay, and and 
are the decay rate and the number density of the heavier particle. In contrast 
to refs. [T2] we assume that does not dominate the total energy density, so that 
the co-moving entropy density remains approximately constant throughout. The 
Boltzmann equation for can then easily be solved analytically, using the fact 
that t oc T~ 2 oc x 2 in the radiation-dominated era. Inserting this solution into the 
equation for n x , and again switching variables to Y x = n x /s, Y$ = n^/s and x, the 
Boltzmann equation for x becomes 



dY v 



(<jv)s 



(Y 2 - Y*J + NrxY^x ) exp [--{x 2 



■In 



dx Rx^ x ~ xm ' ' *v~"/~x- v 2 v- ~»J) > ( 35 ) 

where r = T^/Hx 2 = (r^M P1 / nm 2 ,) -\/90/g* is constant. The zeroth order solution 
of Eq. (j33j) is again obtained by neglecting x annihilation. Using the expansion (jSJ) 
of the annihilation cross section, we have 



dY 



f yi+ — j cxe 2x + NrxY^xo) exp 



dx \ x 
This equation can be integrated, giving 

Yo = fc[^(x e- 2xo -xe- 2x ) + 
For i>i , Y becomes constant 



(36) 



+ 36 



-2x 



— e 



-2x\ 



1 - exp (--(x -x Qi 



L2 



x e 



+ 



36 e 



-2x 



+ NY^xo) + Y x (x ) 



(37) 



(3f 



For sufficiently large Yq the annihilation term in Eq. ()35|) becomes significant. We 
add a correction term to include this effect, as in Eq. (|2()j) . Since the new, non- 
thermal contribution to x production is already fully included in Y , the Boltzmann 
equation for 5 is again given by Eq.(|2T|). Using now Eq.(|3T|) for Y , we can integrate 
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Eq.flSU), giving 



-fc 2 



^-a 3 F 4 (x, x ) + ^a 2 (a + 186)F 4 (x, x ) 

1 3 
+—a(a + 126) (a + 366)F 2 4 (x,x ) + -b(a + 126) 2 F 3 4 ( 
16 8 



a 2 F{(x, x ) + -a(a + 246)F 2 (x, x ) + 36(a + 126)F 2 (x, x ) 



+lo,oo/ 2 c 

[aF 2 °(x, x ) + 66F 3 °(x, x )] 

- N 2 Y*(x o y*of[aG r 2 (x,x ) + 6bG r 3 (x,x )} 

+ 2NY <j> (x )e rx ^ 2 Y , oo f[aG r 2 /2 (x, x ) + 6bG r 3 /2 (x, x )] 



(39) 



- NY^x )e rx o/ 2 fc 



a 2 Gl(x, x ) + -a(a + 246)G2(x, x ) + 36(a + 12b)Gl(x, x ) 



The functions G r „( X, Xo), G-J (x : Xq ) and G^(x, Xq) are denned by 









G r n 


(x,x ) = 


j dt 






J XO 








V X, Xo) 


1 dt 






Jxo 










v X,Xq) = 


1 dt 






Jxo 



-rt 2 /2 
t n 



n 



n 



2,3, 



2,3 



n 



1,2,3. 



(40) 



Explicit expressions for these functions are given in the Appendix, Eqs. (j48)l . Notice 
that the expression in curly brackets {...} in Eq.(|39|) has the same form as in Eq.(|22j). 

Results for this scenario with 6 = are shown in Fig. [HI We choose r = 0.1 so that 
rx\ ~ Xo, which leads to the most difficult situation where thermal and non-thermal 
production occur simultaneously. We see that even for the smaller cross section 
considered, a = 10~ 9 GeV~ 2 (top frames), the simple first-order solution (J20j) soon 
fails, since \S\ exceeds Yq. However, the re-summed ansatz Y\, r of Eq.()25|) describes 
the exact temperature dependence very well for this cross section, both for large 
(top left frame) and moderate (top right) non-thermal \ production. For a = 10 -8 
GeV -2 (bottom frames) we again observe sizable deviations for intermediate values 
of x — Xo- 

In fact, comparison with Fig. |21 shows that non-thermal \ production leads to 
faster growth of and hence to earlier and larger deviation between Y" 1/r and the 
exact solution of the Boltzmann equation (J35|) . However, comparison with the curves 
labeled K^tp, where non-thermal x production is neglected, show that for this rather 
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0.01 0.1 1 10 0.01 0.1 1 10 



x x x x 

(c) (d) 

Figure 8: Evolution of Y x (solid red curves), Y\^ r (dotted blue), \8\ (double-dotted black), the 
prediction for purely thermal x production Y Xttp (short-dashed violet) and Y cq (triple-dotted or- 
ange) as function of x — xq, for Y x (xo = 22) = 0, r = 0.1, N = 1 and 6 = 0. The S— wave cross 
section and the initial 4> density are (a) a = 10~ 9 GeV~ 2 , Y^(xo) = 10~ 10 , (b) a = 10~ 9 GeV~ 2 , 
Y^xq) = 10" u , (c) a = 10~ 8 GcV- 2 , Y^xq) = 10~ 9 and (d) a = 10~ 8 GeV~ 2 , Y^xq) = lO" 10 . 



large cross section and short lifetime, the non-thermal production mechanism does 
not affect the final x re li c density any more. This agrees with the result of Fig. El 
where we saw that for the same values of a and Xq, the relic density is independent 
of the initial value Y x (xq). As before, Yi >r approaches the exact result again for 
x — Xq 3> 1. We therefore conclude that our re-summed ansatz describes scenarios 
with additional non-thermal x production as well as the simpler case with only 
thermal production. 
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5 Summary and Conclusions 



In this paper we investigated the relic abundance of non-relativistic long-lived or 
stable particles x using analytical as well as numerical methods. Our emphasis 
was on scenarios with low re-heat temperature, so that x m ay never have been in 
full thermal equilibrium after the end of inflation. Such scenarios are interesting 
because they lower the predicted relic abundance and therefore open the parameter 
space of particle physics models, allowing combinations of parameters which are 
cosmologically disfavored in the standard high temperature scenario. 

The case of small x annihilation cross section or very low temperature can eas- 
ily be treated analytically, since in this case x annihilation can either be ignored 
completely, leading to our zeroth order solution Y Q of Eq. (JT7j) . or can be treated as 
small perturbation, as in our first order solution Yi of Eq. (|2U|) . Unfortunately this 
approximation breaks down well before x attains full thermal equilibrium. On the 
other hand, we found that the simple trick of "re-summing" the correction due to 
X annihilation, as in Eq. (j25|) . allows to describe the full temperature dependence 
of the x number density as long as x does not reach full equilibrium. We saw in 
Sec. 4 that this remains true even if a non-thermal source of x production is added. 
Our ansatz therefore provides a first analytical description of the "in-between" situ- 
ation, where x annihilation is very significant but not large enough to establish full 
chemical equilibrium with the thermal plasma. 

For yet higher cross sections or temperatures even the re-summed ansatz fails 
to describe the temperature dependence of the x number density at intermediate 
temperatures. However, by replacing the initial scaled inverse temperature Xq with 
the quantity i ,max of Eq. our ansatz succeeds in predicting the final relic density 
about as well as the standard semi-analytical high temperature treatment does, with 
comparable numerical effort. 

In this paper we have used the non-relativistic expansion of the x annihilation 
cross section. This expansion is known to fail in certain cases even for non-relativistic 
WIMPs We expect our methods to be applicable to these situations as well. 
However, a full analytical treatment will be possible only if the product of ther- 
mally averaged cross section and squared x equilibrium number density, expressed 
as function of the scaled inverse temperature x, can be integrated analytically over 
x. 

From the particle physics point of view, the main effect of a low reheat temper- 
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ature is that it allows to reproduce the correct relic density in scenarios with low 
annihilation cross section, e.g. for Bino-like neutralinos and large sfermion masses. 
Conversely, the non-thermal production mechanism studied in Sec. 4 allows to re- 
produce the correct relic density for WIMPs with large annihilation cross section, 
e.g. Wino-like neutralinos 21. As noticed in jT3], the combination of these effects in 
principle allows to completely decouple the WIMP relic density from its annihilation 
cross section. In many studies of expected WIMP detection rates scenarios yielding 
too high a relic density under the standard assumptions were not considered; such 
scenarios typically also lead to low detection rates. Conversely, in scenarios leading 
to too low a thermal WIMP density, which typically predict large detection rates for 
fixed WIMP density, the predicted detection rates were often rescaled by the ratio of 
the predicted to the observed relic density. If one allows lower reheat temperatures 
and/or non-thermal WIMP sources the possible range of signals for WIMP detection 
can therefore be enlarged towards both larger and smaller values. 

In summary, we found analytical or semi-analytical solutions of the Boltzmann 
equation describing the density of non-relativistic relics which are valid for a wide 
range of initial conditions. In particular, they allow a complete description of the 
temperature dependence for small or moderate cross sections, and correctly repro- 
duce the final relic density for all combinations of initial temperature and cross 
section. This should be a powerful tool for exploring the physics of non-relativistic 
relics, especially in scenarios with low reheat temperature. 
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Appendix 

In this Appendix, we give explicit expressions for the functions F™(x, Xo), G r n {x, Xq), 
Gn (x,xq) and G c n [x, xq) which appear in Sees. 3 and 4. These functions are ana- 
lytically expressed in terms of the exponential integral of the first order Ex (x) and 
the error function erfc(x). 
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First we review the exponential integral and the error function. The exponential 
integral of the first order is defined by 

/oo -xt roc -t 

dt — = J dt—. (41) 

We need this function only for x > xq ^> 1. We can then use the asymptotic large x 
expansion, 

e~ x A (-l) n n\ 



(42) 

n=0 

The error function is defined by 

erfc(x) = —= / die"' 2 , (43) 
V 71 " Jx 



with asymptotic large x expansion 



ericW-^f; '- 1 );^- 1 " 1 . (44) 



n=0 

The functions F™(x,xo) are defined by 

rx -rat 

F™(x,x )= dt—. (45) 

Jxq 1 

These integrals can be reduced to the form (|41|). The resulting expressions and 
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corresponding asymptotic expansions, computed from Eq.(jI2J), are: 

F \x,x ) = i(e-^-e- 4 -), 
Ff(x,x ) = Ei(4xo)-Ei(4x) 

Ax \ Ax ) Ax \ Ax J \ Xg 
e -4x e~ 4x 
F 2 4 (x,x ) = 4E!(4x ) + AE 1 {Ax) 

Xq X 
p-4x q~ 4x /p-4x 

+ o 



F 2 °(x,x ) 



4xq Ax 2 \ Xq 
g-4x e -4xo e _4a; e~ ix 

Fi(x,x ) = —-2— + 8E 1 (Ax )- — + 2—-8E 1 (Ax) 

/ p -4x 

F^x,^) = E 1 (2x )-E 1 (2x) 

-2x / 1 \ P ~2x / 1 \ / p-2a:o 



2^o V 7 2i \ 2x J \ Xq 
-2x e _2x 
F%(x,x ) = 2E 1 (2£ ) + 2Ei(2x) 

Xq X 
p-2x p,~ 2x f e~ 2x 

+ o 



2x1 2x 2 



^0 ^ \ x o 

F*(x,x ) = — — + 2E 1 (2x )-^ + ^-2E 1 (2x) 



3 



e ~2x e -2xo e~ 2x e~ 2x 



2x\ Xq 2x 2 X 

-2x ' 

O ' ' 



c 3 

I 1 

Xq X ' 



i^Wo) = 2^-^»- ( 46 ) 
The functions G r n (x, x ) and Gn(x,x ) are defined by 

rx e - rt ' 2 

G r n (x,x ) = / dt——, n = 2,3, 

</x 1 

Gf(i.xo) = / dt e ——, n = 2,3. (47) 

JXQ 1 

Using Eqs. (|43jl and (|44jl . we find the following explicit expressions and corresponding 
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asymptotic expansions: 



-rx 



G r 2 (x,x ) 



a 



.—rx 



G r 3 (x,x ) 



2txq 



2x\ 



2tXq 



~r eric(y/rxo) 
3 



x 



2rx\ 



2rx 3 



3 • o 



2rx 2 



x (rxlY 



o E iK) - 



i - 



rx 



— + - Ei(rx 2 ) 

2x 2 2 v ' 



— * +o 



2rx 4 



^o( ra o) 3 



G2 /2 (a;,xo) 



e r:r o/ 2 /7rr 

\ — erfc 

x V 2 

p -rx 2 /2 



2"° — 




+ v/y erfc 



r:r 



o 



-rx 2 /2 



G r 3 /2 {x,x ) 



-rx 2 /2 

~2xY 

-rxg/2 



rrrf 



4 1 V 2 



3 

r.r 



+ I — - 

- ' ' x (rx^ 3 




e ra2 / 2 r / rx 2N 



1 



"0 \ rx o, 

In the expansion we assume that rx\ ~ Xq, so that the effect of non-thermal \ 
production is comparable to that of thermal production. 
Finally, the functions G„(x, xq) are defined by 

-2t-rt 2 /2 



-rx 1 j2 



1 



4 

rx 



+ O 

- ' ' x1(rx1) 3 



(4f 



xq) 



dt 



n 



1,2,3. 



(49) 



They appear in the "interference terms" in Eq. (j39j) . which are important only if 
thermal and non-thermal contributions to Yq in Eq. ()37|) are comparable in size. 
Since the overall t— dependence of the integrand in Eq. (j49|) is dominated by the 
numerator, we can, to good approximation, evaluate these functions by replacing t 
in the denominator by some appropriate constant x c : 

-2t-rt 2 /2 



G c n (x,x ) 



dt 

2/r 



X : 




7T 

2r 



erfc 



(rx + 2) 



e -2x -rxl/2 

x™(rx + 2) 



(rx + 2) 2 



erfc 



e -2x-rx 2 /2 

x^irx + 2) 



rx 



(rx + 2) 2 



-2x -rx 2 /2 



xr\rxlf 



(50) 
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In our calculations in Sec. 5 we set x c = xo; this over-estimates G c n by a few %, with 
negligible error in Y" l r . 
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